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ABSTRACT 



We derive and investigate the dispersion relation for accretion disks with retarded or advanced heating. We follow the 
«-prescription but allow for a time offset r between heating and pressure perturbations, as well as for a diminished 
response of heating to pressure variations. We study in detail solutions of the dispersion relation for disks with radiation- 
pressure fraction 1 — /3. For r < (delayed heating) the number and sign of real solutions for the growth rate depend 
on the values of the time lag and the ratio of heating response to pressure perturbations, ^. If the delay is larger than a 
critical value (e.g., if Qt < —125 for a = 0.1, P = and ^ = 1) two real solutions exist, which are both negative. These 
results imply that retarded heating may stabilize radiation-pressure dominated accretion disks. 

Key words, black holes physics — accretion disks — stability — time-delay 



1. Introduction 



The IShakura fc SunvaevI (|1973l ) a -viscosity accretion disk 
model has been extremely successful in describing vari- 
■ ous astronomical objects and systems. The only excep- 
' tion is its application to systems accreting at high rates. 
At rates where pressure is dominated by that of radiation 
and opacity by electron scattering, t he a disk is thermall 
and viscously (secularly) unstable (iLightman fc Eardlc 
[T97I IShakura fc SunvaevI . [19761 : iShibazaki fcHoshj llQT^ 
In the case of accretion onto black holes this regime cor- 
responds to luminosities in excess of > 0.01 iEdd, where 
-^Edd = GMrUp c/aT, Trip is the proton mass and ut the 
Thompson cross-section. However, black-hole X-ray sources 
cross this limit upwards to maximum luminosity and down- 
wards to minimum luminosity showing no dramatic symp- 
toms at all (but they enter the so-call ed hard/low state, 
see e.g., iMcClintock fc Remillardl . I2006D . and certai nly not 
the behavior antic i pated by models (e.g., |Lasota fc Pelatl . 
I1991I : iTaam fc Linl . [1983) ■ Observations suggest that disks 
in black-ho le transient system s are stable up to at least 
~ 0.5LEdd ([Done et al.l . l2004[ ). 



Models of radiation-pressure dominated disks are unsta- 
ble only when the viscous stress is proportional to the to- 
tal pressure {a being the proportionality constant). Models 
with ad hoc viscosity prescriptions have been studied but 
their relation to reality remains unknown. Simulations 
of the magneto-rotationa l instability in radiation-pressure 
dominated disks (see e.g. . [Turner! . 1200 J) showed that stress 
was approximately proportional to the total pressure, but 
they exhibited no sign of instability. 

Recently, iHirose et al.l (|2009l ) showed that although the 
linear correlation of vertically integrated stress and pres- 
sure is roughly satisfied in shearing-box MHD simulations 
of radiation-pressure dominated disks, these quantities are 
shifted in time — pressure responds to stress variations after 
10—20 dynamical times. Using these simulation results as 
a guideline, we perform an analytical, perturbative study of 
the stability of disks with a modified viscosity prescription 
that allows for such a time lag between stress and pressure. 

The theory of dela yed oscillators has already been exten- 



Minorskv[ . ll944t[Cooke fc Grossmanl . 



sively developed (e.g., , . _ ^ _ ,, 

Il98a iBellman fc Cookd . Il96.tl . It predicts that an oscilla- 
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tor may be easily stabilized (or destabilized) if only the 
system parameters are chosen properly. We show that this 
is the case also for accretion disks. Similar conclusions were 
recently obtained in anot her analytic study parallel to this 
work bv lLin et all (|201lD . 

This work is devoted to the mathematical part of this 
project. In Sect. 12.11 we discuss our choice of the modified 
viscosity prescription. In Sect. 12.2] we derive the dispersion 
relation. In Sect. [3] we present a detailed discussion of the 
long wavelength limit. In Sect.H]we discuss solutions of the 
dispersion relation. Finally, in Sect. [5] we summarize our 
results. Detailed discussion of their physical implications 
will be given in a separate paper. 

2. Perturbative analysis of disk stability 

We base our study on a linear pe rturba t ive an alysis fol- 
lowing the approach pioneered by iPiranI (|1978[ ) . The un- 
perturbed disk is assumed to be steady, so that time de- 
pendence can only be found in the perturbations. In the 
following, it is understood that the characteristic length- 
scales of radial variation of unperturbed variables are ^ r. 
We consider small axisymmetric perturbations with radial 
wavenumber fc, assuming that their wavelength A = 27r/fc 
satisfies the relation 

i?<A<r, (1) 

where H is the half-thickness of the disk. (The assumption 
of geometrical thinness ceases to be valid for disks with 
L > few X 0.01 LEdd, so that strictly speaking our results 
apply only to the low luminosity region of the regime where 
the thermal-viscous instability would appear according to 
the standard model). Since the radial flow is slow in com- 
parison to the azimuthal m otion, we make the usual as- 
sumption (|Kato et al.l . |2008[ ) in calculations of instability 
that Vro = 0, where the index will denote unperturbed 
quantities. We define the following dimensionless variables 
corresponding to the Eulerian perturbations of vertically in- 
tegrated pressure (Pi), radial velocity (vri), surface density 
(Si), disk thickness {Hi), and vertically integrated viscous 
stress (Ti), 

ilr So Po Ho To 

where n = ^jGMjr^ is the Keplerian rotational frequency. 
We assume that all of these quantitites represent complex 
waveforms e"^*^*'^'', e.g., vj = tue"^*^*'^'' and 

P{v,t) = Po{r) + Pi(r, t) = Po(r) • (1 + ^e""*"'^'^^), (3) 

where n stands for the dimensionless frequency, and w is 
a constant and uniform dimensionless amplitude. Negative 
values of the real part of the dimensionless frequency, K(n) , 
correspond to damped (stable) perturbations while positive 
values to exponentially growing (unstable) ones. The imag- 
inary part of n determines the frequency of the correspond- 
ing oscillations. 

2. 1 . Viscosity prescription 

In the standard approach , based on the a-prescription 
(jShakura fc SunvaevL [l973) , one assumes that the r(p com- 
ponent of the stress tensor, T^^, is proportional to pressure: 

Tr^ir,t) ^ -aP{r,t), (4) 



where Tr^^ and P are given in terms of the unperturbed and 
perturbed quantities as 

Trv =ro(l + ee"''*"*'^'''), (5) 

P = Po(l + n7e""*~*'=''). (6) 

In this work we assume a modified prescription for vis- 
cosity. Instead of assuming that a is constant both in time 
and radius we write 

a{r,t) =ao{l + a) = ao [l + (^e""'^" - tu)e"^"-'^'-] , (7) 

with 

To = -aoPo. (8) 

As 

-a{r,t + T)P{r,t + r) = 

xPo(l-f we""(*+^)-''='') 
« -aoPo(l + 0e""*-''^'-), (9) 

taking into account Eqs. ([5|) — (O we have, through first 
order in the perturbation, 

Tr^{r,t) ^ -a{r,t + T)P{r,t + T). (10) 

Hence, our modified viscosity prescription, Eq. (I7|), corre- 
sponds to a stress which is still proportional to aP but 
is advanced with respect to this quantity by the time r, 
i.e.. the value of stress at time t is proportional to aP as 
measured at time t + t. In particular, a negative value of 
T corresponds to a delayed stress: a perturbation of stress 
follows the perturbation in pressure after a delay of — r. 

We also introduce the heating to pressure response fac- 
tor, defined as the ratio of the dimensionless amplitudes of 
stress and pressure perturbations. 



We assume ^ to be real. 

2.2. Derivation of the dispersion relation 

In the following four subsections we derive perturbed forms 
of hydrostatic balance, mass conservation, angular momen- 
tum balance and energy equation for geometrically thin, 
axisymmetric accretion disks. 

2.2.1. Hydrostatic balance 

The balance of vertical forces, after integrating along the 
vertical coordinate, takes the following form, 

n^H^ = ^. (12) 

Writing this form of the vertical force balance we assume 
the disk to be in hydrostatic equilibrium, since the thermal 
and secular timescales are much longer than the dynamical 
timescale. 

Perturbing Eq. (jl2p with small-amplitude axisymmetric 
perturbations, and assuming that the azimuthal component 
of the velocity undergoes no change, we get 

^'H^u{l + W = ^%^. (13) 
So(l + cr) 
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Using the unperturbed Eq. ([T^ . we get 

{l + h)^{l + a) = l + w, 
i.e., through Hnear order 

l + 2h + a^l + w, 

and finally 

2h = m — a. 
2.2.2. Mass conservation 



2.2.4. Energy equation 

We start f rom the second l aw of thermodynamics in the 
form fe.g.- lKato et al.l . l2008D 



dE_ 
'dt 



P- 



ot ror or 



-VrP- 



d\nH 



dr 



)+ -Q- 



(14) or simply. 
dE d 



dH 



rdr 



{rVrH) 



+ 

vis 



The vertically mtegrated form of the contmuity equation , n ■ ii i- n ■ ^ ^ j -n 

T ■,, J ^ where E IS the vertically integrated specmc energy, 

can be written as j a t- oj^ 



9S 
'dt 

Its perturbed form is 



+ ^(rStv) =0. 
ror 



— Eo(l + a) + ^[rSo(l + ij)rnu] = 0. 
ot ror 



(15) 



(16) 



Neglecting terms of the second order (ua) and using 
Eq. (fT5)l in its unperturbed form, dYlo/dt = 0, we obtain 



din (r^Sof^) 
91nr 



ikr 



= 0. 



(17) 



The logarithmic derivative term may be neglected on the 
strength of the assumption stated in Eq. ([1]). Thus, we ob- 
tain the final relation 



na = iKru. 



2.2.3. Angular momentum conservation 

The angular momentum conservation law is 

d 



Evr--^{r^n) 
r dr 



r^d'i 



-{r'Tr^). 



(19) 



For the first order perturbations we obtain 

ll,^ru~{r^n) = -\^{r^T,i 
r dr dr ^ 



According to Eq. ([T|) the derivative on the right hand side 
is dominated by the dO/dr term. Using Eq. ([5]), and intro- 
ducing the speed of sound squared = Pq/^Oj ^-nd the 
vertical epicyclic frequency 



2r!d(r2r!) 



dr 



we obtain 



and finally. 



2172 



0, 



u = 2ikrao^ I 



V7 



(21) 



(22) 



(23) 



Incidentally, in our study k, ~ Q and Eq. (jl2p reads 
Cs/{rn) = Ho/r, so Cs/irn) = Ho/r. 



E 



3(1-/?) + 



7-1 



P = AP, 



(24) 



Qrad^ 

(25) 
(26) 



P is the gas to total pressure ratio, 7 is the ratio of specific 
heats, Qvis ^^'^ Qrad ^^'^ viscous heating and radiative 
cooling rates per unit area, respectively. 

Let us differentiate both sides of Eq. (|25p with respect to 
time. Assuming that Q^^g is a function of P and a (following 
the a-prescription) while Q^^^ is a function of P and E 
(as is the case for radiative cooling in the optically thick 
regime), we obtain, 



dt 
dt 



dP 
dP 



dP 

dt 



da 



dt \ as 



9S 
'dt 



(28) 



(18) Introducing the perturbations and differentiating we get. 



1 dQL ^ (dOL 
nfl dt V 9P 

±dQrad ^ (dOrad 

nn dt V dP 



PqW - 



Pom 



E.Q 



da 
\ dT. 



aod (29) 



SofT (30) 



p. a 



Finally, we have, 

^ {Qt^s - Qrad) = Po i^G^ + + ««„) , (31) 

(20) where 



Gcr 
Ga 



1 / fdQt:. 



fdQ- 



n U dP 

1 (dQ^ao 



\ dP 



VLcl \ dT. 



ao 



nPo V da 



(32) 
(33) 
(34) 



Following the g-p rescription, these quantities simplify to 
(|KatoliaD,[200i), 



Gtu = 

Gcr = 
Ga ~ 



3ao (2 - 5/3) 
~ ■ (4 - 3/3) ' 
3ao (2 -t- 3/3) 
2 ' (4 - 3/3) ' 
3ao 



(35) 
(36) 
(37) 
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T he perturbation of is given by the following expres- 
sion (|Chandrasekaii [19581) ■ 



Po 



l-/3o 
4-3/3o 



(4ct - 3tt7 - h). 



(38) 



Hence, like all the other time derivatives, the time deriva- 
tive of A is first order in the perturbation: 



dA 
'dt 



d_ 

dt 



3(1-/3) + 



4 - 37 

' 7- 1 at ■ 



(39) 



7-1. 

The time derivative of the first term on the left hand side 
of Eq. ((25|) is given by. 



d^AP _ ^d^A 



dA dP 



d^P 



a<2 df^ ^ df^ ' ^ dt dt ' ^ dt^ ' ^^^^ 

where Eq. has been taken into account. The derivatives 
of A and P are both first order, so the middle term is second 
order and may be neglected. The other terms are 

d'^P 



A 



and 



df^ 
d^A 



= P 



4-37 



dt^ 7-1 dt^ 
The second time derivative of /? = (3o{l + (3) equals, 
d^P 



2o2/ 



df^ 



so finally, 



d^E 



= nfl P() nAw 



4-37 



7 



1 



(41) 
(42) 

(43) 
(44) 



Perturbing the second term on the left hand side of 
Eq. keeping in mind the assumption of Eq. ([T]), we 

get through first order 



-^{r^nuE) 
rdr 



-ihrflAPou. 
The time derivative of this expression is 

-ikrnil^PoAu . 



(45) 



(46) 



The time derivatives of the terms involving H in 
Eq. (|25p are straightforward to compute. The first order 
terms are 



P 
H 



dH 
'dt 



d 
rdr 



(rVrH) 



PonQh 



Po d 



Ha rdr 
P^nQh — ikPorflu 



(r^nuHo) 



(47) 



and their time derivative is 
d 



dt 



PdH p d , 



^nV?PQ{nh-ikru). (48) 



Collecting Eqs. ([311) , (El) , dS]) and (gS]) , we obtain the final 
form of the perturbed energy equation. 



nAuj+n 



4-37 
7-1 



l3QP+nh-ikru{A+l) = wG^+aGa+aGa- 

(49) 



2.2.5. System of perturbed equations and the dispersion 
relation 

Eqs. ((131), (O, (221) and (03) form a system of four algebraic 
equations. 

The perturbation of a, defined in Eq. ([71), may be ex- 
pressed in terms of w, and ^ of Eq. (jlip : 



a 



(50) 



Taking this relation into account we obtain four coupled 
homogeneous algebraic equations in the following form, 



2h = Tu — a , 
na = ikru, 

u = 2ikrao 



(51) 



AnvD 4 
= wG-n 



i^^PoP + nh- ikru{A + 1) 



where /? is given by Eq. (|38p . This system has nontrivial 
solutions for the four variables /i, -kj, ct, and u, iff 

n'Ci + n [^FC2 + G'„(l - ^e"""") - G^] + ^FG^ = 0, 

(52) 

where 



F 
Ci 

B 



2ao 



A 
A 



^ ^(37-4)(l-/3) 



(53) 
(54) 
(55) 
(56) 



(7-l)(4-3/3)^ 

and A is given in Eq. ([^5)) . For (3 = the coefficients of 
Eq. ([5^ do not depend on 7. For ^ = 1 and fir = the 
condition given in Eq. (j52p simplifies to the stan dard disper- 
sion r elation for perturbed accretion disks (^^e.g.. lKato et al.l . 
I2OO8D . Note that < ^ < 1, and that for 7 > 4/3, A and 
B satisfy A > 0, and B > 0. Thus, both Ci and G2 are 
positive for any value of (3. In the following, we will take 
7 = 5/3 whenever a specific value is required for numerical 
results. We also specialize to F = 2aQ{kHo)^ , in accordance 
with the comment following Eq. (|^5)) . 

3. Long-wave limit 

Let us first consider the solutions of the dispersion relation 
Eq. in the long wavelength limit, i.e., kH — > 0, which is 
very useful in classifying solutions of arbitrary wavelength. 
In the limit kH — !> we neglect terms proportional to F 
and obtain 

n [nCi + (1 - Ce-""")G„ - G„] = 0, (57) 

with a trivial solution n = 0. Dividing by n 7^ and using 
Eqs. ([55)) - ([57)) we get for the remaining solutions 

3 



where 



"Gi + -ao (1 - ee-"''" - /) = 



G^ ^ 2-5(3 
■' ~ G„ 4-3;3' 



(58) 
(59) 
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1.0 - 




Fig. 1. The shaded area denotes the region in the (/?, ^) 
plane for which the real root of Eq. ([55]) is negative for all 
Qt > 0. 



3.1. The case of no delay, fir = 



Solving Eq. (|58p is easy for Or = 0. In this case n is real 
and is given by, 



no 



3ao(C + /-l) 
2Ci 



(60) 



As Ci is positive for all values of /3, the sign of n is deter- 
mined by the sign of ^ + / — 1. For ^ = 1, no = Gt^/Ci 
and it is positive for /? < 2/5 and negative for /3 > 2/5. The 
latter inequality is the standard condition for disk stability. 

In general, we get the following criterion for the negative 
sign of n, 

^<-/^^- 

The shaded area on Fig. [1] presents the region in the (/3, 
plane for which n is negative. Note that it suffices to de- 
crease the amplitude of stress variations by a factor of two 
= 0.5) to stabilize the disk for all /3. Despite the fact 
that the condition (|FT|) has been derived assuming VIt = 0, 
it remains satisfied for the negative roots of Eq. ([55)) for all 
VIt > 0, as will be shown in Section [ 



3.2. Arbitrary fir 

Let us now consider the general case of fir ^ 0. 
Equation (j58p is no longer trivial as it involves an exponen- 
tial function of n. An infinite number of complex solutions 
is expected as the exponential function e""^^"^ involves pe- 
riodic trigonometric functions whenever the imaginary part 
of n is nonzero, 3(7i) ^ 0. All imaginary solutions are conju- 
gate, i.e., if X satisfies Eq. ([55)) then x* is also a solution. As 
will be shown in the following section, no more than two 
real solutions may exist. To find the roots of a nonlinear 
complex equation such as Eq. ()58p one has to use numer- 
ical m ethods. We used the MINPACK routines ([More et all . 
|1984[ ). First, the locations of minima of the absolute value 
of the left hand side of Eq. ([55|) were roughly estimated. 
The values so obtained served as starting points for the 
nonlinear solver. 

In Fig. [2] we plot color-coded absolute values of the left 
hand side of Eq. ([55)) for /3 = 0, ^ = 1, a = 0.1 and three 
values of fir = -50 (left), = (middle), f7r = 50 



0.025 - 



-0.025 







Fig. 3. Real part of solutions of Eq. {[55|) for /3 = 0, ^ = 1 
and a = 0.1 as a function of r. Tire solid lines denote real 
solutions, the dashed line shows complex solutions linking 
the real branches, while the dotted lines present periodic 
complex solutions (only first 20 solutions are drawn) . Stars 
denote points where the real branches merge into the com- 
plex one (and its conjugate). 



(right panel). Only 3(n) > regions are shown. The darker 
the color, the smaller the value. Red crosses denote real 
solutions while red squares show locations of solutions with 
non-zero imaginary part. For each of the values 51r = 
and fir = 50 a single real solution exists, and it satisfies 
3fJ(?i) > 0. For fir ^ there is an infinite number of complex 
solutions. The sign of their real part (with the exception of 
the first imaginary root when fir < 0) is in general opposite 
to the sign of fir. 

Fig. 121 presents the real part of solutions of Eq. ([55)) 
in the (r2r,5R(n)) plane for the chosen values of a, /? and 
^. Solid lines present real solutions (crosses in Fig. 
Dotted lines correspond to ordinary complex roots (squares 
in Fig. (2) — only the first 20 are plotted. The dashed line 
connecting two real branches is a special class of complex 
solutions. 

Taking fir = we recover the real solution giveir by 
Eq. ([50)) . For small but negative fir there are two real so- 
lutions, one of which satisfies lim^^^Q- 3fi(n) = -|-oo. These 
two real solutions converge to each other and merge into 
complex conjugate solutions (empty star in Fig. [5]) at a 
critical value of firi « — 11 . These roots split into two 
real solutions, but this time negative, at 17x2 ~ —125 (solid 
star). The upper branch of these two approaches 5R(n) = 0~ 
for fir — — oo. For fir > there is only one real solution, 
which is positive and approaches 3?(7i) = 0+ for fir +oo. 

In addition to these solutions, an infinite number of pe- 
riodic complex solutions exists (dotted lines). For most of 
the range of fir presented in Fig. |5|they appear in the sec- 
ond and fourth quadrant of (fir, K(n)) plane for negative 
and positive fir, respectively. For |fir| > 150, however, the 
first complex root, and subsequently the others too (but at 
much larger absolute values of fir), crosses the 3fi(n) = 
axis. 
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Fig. 2. Maps presenting absolute values of the left hand side of Eq. [SHI (dispersion relation in the limit of long waves) 
on the complex plane for /3 = 0, ^ = 1, a = 0.1 and three values of fir = —50 (left), (middle) and 50 (right panel). 
Crosses and rectangles denote locations of the real and complex solutions, respectively. 



3.3. Real roots 

Let us examine the real branches in detail. Eq. (j58p may be 
rewritten in the following form, 



instance, under which Eq. (j58p has only one real solution. 
After some algebra we get. 



ce 



where 
a 
b 



Ci>0, 

^ao(l-/)>0, 
3 

-ao^ > 0. 



(62) 

(63) 
(64) 

(65) 



b_ 

ce 



1-/ 



we 



where 



^^,_ 3ao(l-/) 
a 2Ci 



(67) 



(68) 



Convexity of the exponential function on the right hand 
side of Eq. (|62|) implies that there are no more than two 
real numbers n satisfying this equation. For Qt > the 
exponential function is decreasing with n and therefore has 
exactly one intersection with an + b, and only for n > (or 
n = 0, or n < 0) as long as b/c < 1 (or 6/c = 1, or 6/c > 1, 
respectively). This yields the same criterion for stability as 
in the fir = case, inequality (j6ip . In this sense, advanced 
heating [T'r0(i) = ~aP{t -\- r) with r > 0] does not alter 
the viscous and thermal stability properties of the disk. 

For fir < the exponential function is, on the contrary, 
increasing and may not intersect the linear one at all, may 
have a single intersection point or may cross it twice. All 
three cases occur and are clearly visible in Fig. [31 the sin- 
gle solutions are marked by stars and the complex branch 
connecting these real solutions corresponds to the no inter- 
section case. 

The condition for a single real solution of Eq. ([55)) for 
fir < may be obtained by matching the gradient of the 
linear function with the derivative of the right hand side of 
Eq. ([SH), 



Solutions of Eg. ([571) are giv en by the multivalued Lambert 
W function ([Wrightl . |1959( ) of its left hand side. Real so- 
lutions exist only for values of the left hand side of the 
equation greater or equal than — 1/e. This condition corre- 
sponds to, 

2(1 + /?) 



^ > 



A-ij3 



(69) 



As we^ has a single minimum, of value —1/e at w = — 1, 
Eq. ([67| has two solutions for all —1 < w < 0, i.e., when 
the inequality in ([59| is sharp. As long as condition (|69p is 
satisfied, the values of fir for which Eq. ([55)) has only one 
real solution are therefore given by 



-/i.2(e,/3) 



1 

"0 ' 



(70) 



dn 



(66) 



Eqs. ()62p and ()66p form a set of two equations correspond- 
ing to the linear function being tangent to the exponential 
one at a given n and providing the condition, on ilr for 



where /i,2(Cj/?) double- valued. We take /i < /2, i.e., 
57t2 < VIti < 0. In Fig. [5] we plot values of /i,2(^,/3) as a 
function of (3 for various values of ^. For the standard case 
of ^ = 1 and radiation pressure dominated disks (/? = 0) we 
have /i(l,0) « 1.08 and /2(1,0) « 12.50. Thus, for, e.g., 
a = 0.1, Eq. ([58p has only one real solution for ftri w —10.8 
and « —125.0 (corresponding to the stars in Fig. [3]). 
A moment's reflection leads to the conclusion that there 
are two positive roots of Eq ([55)) for fJri < fir < 0, two 
negative ones for fir < il,2T, and no real roots for < 
Ut < Uti. When condition ([CT)) is satisfied, instead of ([S^ . 
there are two solutions for n, one positive and one negative. 
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Fig. 4. /i.2(^, /3) (Eq. 170)) dependence on fi for various ^. 
In this figure, we take 7 = 5/3. 



3.4. Limits of infinite time offset (VIt — ±00 j 

Let us now investigate the limits of real solutions of 
Eq. ([55)1 . Assuming that n is finite, Eq. ([5^ simplifies in 
the limit of nfir — > +00 to. 



3ao 2(1 + 
2Ci 4-3/3 



< 0. 



(71) 



Therefore, this value is valid only for Q.T — >■ — cxd. For the 
standard choice ofa = 0.1, /3 = and C = 1 one gets, 



= -0.021, 



(72) 



which corresponds to the limit of the lower real branch in 
the third quadrant of Fig. [31 

To find limits of the other two branches let us assume 
that nU,T remains finite for \^t\ — > 00. This assumption 
implies 71 — !■ and Eq. (j62p takes in this limit the following 
form, 

fe = ce-"^^^. (73) 
Therefore, n has to fulfill the relation. 



1 1 ^ 
n = -—log-. 

\It c 



The sign of the logarithm depends on ^ and /?: 

2(1 + /3) 



log6/c < for ^ > 
log&/c > for ^ < 



4-3/3 ^ 
2(1 + ^) 
4-3/3 ■ 



(74) 

(75) 
(76) 



Thus, for ^ and /3 satisfying ([75[) n approaches 0"*" for Vtr 
+00 and 0~ for fir — > —00. If condition (|76p is satisfied we 
have lima^^+oo = 0" and hm^T^-oo n = 0+. 



0.05 



-0^^ = 0.02 
-a„ = 0.1 
-a„ = 0.2 



-100 
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£2t 
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100 



Fig. 5. Real part of solutions of Eq. ([55)1 for /3 = and 
^ = 1 as a function of t. Solutions for three values of a are 
presented with different colors (increasing values of ag from 
left to right). Solid lines are for real solutions while dotted 
and dashed lines show the real part of complex solutions. 
Only the first periodic exponential root is plotted. 

curves have qualitatively similar shapes — the sign of solu- 
tions in a given region does not depend on ap (Table [1]) . 
However, the values of firi and riT2, which limit the re- 
gions with double real solutions for fir < 0, are sensitive 
to the value of a^. Eq. ([70[l predicts that their values are 
inversely proportional to ao and therefore, e.g., the region 
with two negative roots of Eq. ([55[) extends to larger values 
of VIt (closer to Or = 0) for high values of a. The limit for 
the lower real branch (n_) also depends on ao according to 
Eq. ([7T|) — the higher the value of a^, the lower the value of 

The impact of /3 on solutions of the dispersion relation 
for a = 0.1 and ^ = 1 is presented in Fig. [Bj For f = 1 the 
crucial inequality ([CT|) corresponds to /3 > 2/5. For values of 
/3 smaller than this critical value solutions exhibit qualita- 
tively the same behavior as discussed previously. The values 
of riTi_2 and n_ depend on (3 according to Eqs. ([70[l and 
([71]) . respectively. Once /3 exceeds 2/5 the character of the 
solution changes. For VIt > the root is negative and ap- 
proaches 0~ with Q,T — > +00. For any negative value of fir 
there are two real solutions with opposite signs — the com- 
plex conjugate branch connecting the real solutions does 
not appear and ilTi^2 are not defined. 

Very similar behavior is shown in Fig. [7) which presents 
the impact of ^ for a = 0.1 and /3 = 0. For this case inequal- 
ity ([CT|) is not satisfied for ^ > 1/2. The solutions change 
their nature once ^ goes below this value, similarly to so- 
lutions with /3 < 2/5 discussed in the previous paragraph. 
In accordance with Eq. (j7ip , n_ , the limit of the lower real 
branch at VLt ^> — 00, does not depend on ^. 

In Table [T] we summarize general features of the real 
solutions of Eq. ([55[) which have been derived in this section. 



3.5. Parameter study 

In Fig. [S] we show the roots of Eq. ([55[) for /3 = 0, ^ = 1 
and three values of ap = 0.02 (green), ao = 0-1 (blue) 
and ao = 0.2 (red line). The second case corresponds to 
Fig. [21 For clarity, from among the infinite number of com- 
plex periodic solutions (dotted lines) , only the first one (in 
the sense of the smallest modulus value) is plotted. All the 



4. Solutions of the dispersion relation 

4.1. Sliort-wave limit 

Let us now consider the short-wave limit {kH ^ 00) of 
Eq. ([5^ . Strictly speaking, this limit violates the first of 
the assumptions of Eq. ([1]). If we assume that n is finite 
then the terms with F dominate. 
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Qt 




fir < 


Or = 


Or > Or - 






2(l+/3) 
i-3/3 


ni — > 


n- <0 
0" 


two negative roots for Qt < Qt2 
complex conjugate roots for f2r2 < Sir < Slri 
two positive roots for ilri < Or 


n = 710 > 


ri > 71 - 






2(l+/3) 
4-3(5 


ni — 
n2 


n_ < 
^0+ 


two roots with opposite signs 


71 = no < 


71 < 71 - 


^ 0" 



Table 1. The long wavelength limit: general characteristic of real solutions of Eq. ((58|) . The quantities uq, 51ri^2 and 
?i_ are given by Eqs. ([50]) . ([70]) and ([7T|l . respectively. 
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Fig. 6. Real part of solutions of Eq. ([55]) for a = 0.1 and 
^ 1 versus r. Solutions for different values of /3 are pre- 
sented (for the solid lines in the fir > region, /3 increases 
from top to bottom) . Solid lines are for real solutions while 
dotted and dashed lines show the real part of complex so- 
lutions. Only the first periodic exponential root is plotted. 
In this figure, we take 7 = 5/3. 
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Fig. 7. Same as Fig. [5] but for a = 0.1, (3 = and varying 
^. For the solid lines in the fir > region, ^ decreases from 
top to bottom. 



n^FC2 + (FGa 



0. 



(77) 



Despite the fact that n in general is complex, in this case 
it must be real to satisfy this real equation. We obtain. 



(78) 



Both G„ = (3ao/2)(2 -f 3^)/(4 - 3^) and Ci are positive. 
Thus, this value of 71 is negative and satisfies the dispersion 
relation in the short-wave limit for all values of parameters 
^ and fir. 

Let us now assume that the absolute value of n is large 
(|7T,| > 1 and \n\ > Ga/C2). In this (and kH > 1) limit 
Eq. ((5^ reduces to 



nCi + iFC2 - Ga^e""^^" - 0. 
Assuming in addition that rifir > we get 
nCi + £,FC2 = 0, 

which is satisfied for 

iFC2 



(79) 



(80) 



(81) 



Hence, this limit is valid for fir < 0. 

If, on the contrary, we assume nflT < then Eq. ([7^ 
reduces to, 



(FC2 - Go.^e 
The solution of which is, 



= 0. 



71 



FC2 

Ga 



Thus, we have two additional limits, 

n — > —00 for Qt > 0, 



n 



-00 for Or < 0. 



(82) 



(83) 



(84) 
(85) 



Summing up, there is a common limit n ~ —Ga/C2 for 
all values of 17t. In addition, solutions with fir > have 
another branch approaching —00 while for negative fir two 
branches are expected approaching both +00 and —00. 

4.2. Arbitrary wavelength 

In this section we consider solutions of the dispersion rela- 
tion for an arbitrary value of the wavelength A = 27r/fc. 

We start with discussing the standard case with no time 
lag (fir = 0) and ^ = 1. Fig. [5] presents the real part of 
solutions of Eq. ([5^ obtained assuming a = 0.1 and various 
values of /3. The limit of l/{kH) 00 corresponds to the 
long- wave limit discussed in Sect. [31 

For /3 < 2/5 there are two positive real solutions: one of 
them approaches zero for long wavelengths (the trivial solu- 
tion of Eq. [S7]) while the other corresponds to no given in 
Eq. According to the classical theory of disk instabili- 
ties they are related to the secular and thermal instabilities, 
respectively. The branches corresponding to the secular and 
thermal modes approach each other with decreasing wave- 
length and merge into complex conjugate solutions (dashed 
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Fig. 8. Real part of solutions of Eq. ([5^ as a function of 
the wave length l/{kH) for ao = 0.1, ^ = 1, fir = and 
various values of /?. The line convention is the same as in 
previous figures. In this figure, we take 7 = 5/3. 
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-0.01 
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Fig. 9. Real part of solutions of Eq. ([5^ as a function of 
the wave length \/{kH) for ao = 0.1, /3 = 0, Ht = and 
various values of ^. The line convention is the same as in 
previous figures. 



lines). Two real and negative branches appear again for 
short wavelengths. The \/{kH) — > limit of the upper one 
corresponds to Eq. (|78p . The lower branch approaches —go 
according to Eq. ([5T|) . The same picture holds for /3 > 2/5, 
except that now all real solutions are negative, at large 
wavelengths one of the roots corresponds to the trivial one 
and approaches 0~ in the limit of 1/ikH) — > cxd, while the 
other one approaches ng < 0. 

For /3 = 2/5 there are two negative real roots at short 
wavelengths, F > 4G^Ci/C2— in the hmit l/{kH) 
one root corresponds to Eq. ([75)) . — G0-/C2 = —0.038 for 
7 = 5/3, while the other tends to -00. At F = 4GctCi/C2 
these merge into complex conjugate solutions tending to the 
origin of the complex plane (n = 0) in the long-wave limit 
{F — 0), while never becoming real for F < AGcrCi/C2 — 
the real part of these conjugate solutions is negative. 

The critical inequality (pT|l relates (3 and ^. Thus, for 
fixed /? similar effects to the ones discussed above may be 
obtained by varying ^. In Fig. [S] we plot roots of the dis- 
persion relation, Eq. ([5^ . for q;o = 0.1, /3 = 0, J7r = and 
a few values of ^. For radiation pressure dominated disks 
(/3 ~ 0) the critical value is ^ = 1/2. For ^ higher than this 
value, the solutions exhibit similar behavior to those dis- 
cussed above for /3 < 2/ 5 — two positive roots for long waves 
(one approaching 0"*"), a common complex conjugate branch 
and two negative real solutions for short waves. The latter 
approach —00 and the limit defined in Eq. (j78p . which does 
not depend on ^. For ^ < 1/2 both solutions are negative for 
all wavelengths and approach the same limits for the short- 
est waves as before. The complex conjugate branch does 
not appear at all for the lowest presented value of ^ = 0.01. 

Fig. (TU] presents similar plots for ao = 0.1, /3 = 0, ^ = 1 
and a few values of the time delay VLt. The long- wave limit 
corresponds to the solutions presented in Fig. [31 The bot- 
tom plot zooms in the shaded region in the top panel. For 
J7r = (black curves) we recover one of the standard cases 
presented in the previous plots. Positive values of fir (e.g.. 
magenta curves) result in two positive roots (one approach- 
ing 0+) for long waves. They merge into the complex con- 
jugate branch with decreasing wavelength and again split 
into two negative real solutions, similarly to some of the 
cases discussed above. 



The behavior of solutions obtained with time delays 
{Q,T < 0) is more complicated. The number of solutions 
in the long-wave limit depends on the relation of the time 
delay to the quantities firi and ^lT2 (for the case presented 
in Fig. [TO] these are —11 and —125, respectively, Eq. \70\). 
For negative Ht > flri we expect in total three real roots: 
two positive, and one equal to zero in the long wavelength 
limit (the trivial one). The red curve in Fig.[TO](ilT = —10) 
corresponds to this case — there are three positive roots, two 
tending to n « 0.04 and w 0.12, and one approaching zero. 

When ilr < $1x2 (e.g., ilr = —150, green lines in 
Fig. [TO)) there are two negative real roots of Eq. ([55]) (com- 
pare Fig. [3]) clearly visible in the bottom panel of Fig. [TOl 
The dotted green line approaching 5R(n) « 0.002 in the 
long-wave limit corresponds to the first complex periodic 
solution (dotted lines in the second quadrant of Fig. [3]). The 
third, positive, solution becomes the trivial one, 3fJ(n) — > 
for kH 0. 

The branches corresponding to the trivial solution of 
the long-wave limit leave the 5R(n) = axis and reach posi- 
tive values with decreasing wavelength. For Qt > ilri they 
merge with the other positive real branch corresponding 
to the smaller of the real roots and transform into com- 
plex conjugate branches (red dashed line). The larger real 
and positive root diverges with decreasing wavelength, ap- 
proaching -l-oo according to Eq. ([55]) . 

Once f2r2 < fir < ftri no real solution of Eq. ([55)1 
exists. The dotted blue line (corresponding to ftr = —12) 
reflects the complex conjugate solution. At large (but finite) 
wavelengths, the single positive solution of Eq. ([5^ corre- 
sponds to the trivial solution ri = of the long wavelength 
limit, Eq. ()57p . However two additional real and positive 
solutions appear at some moderate range of wavelengths. 
These new branches behave similarly to the case previously 
discussed: one of them merges with the trivial branch, the 
other diverges at zero wavelength. 

For Ht < (e.g., green line) the trivial branch di- 
verges on its own, there is only one positive root for all 
wavelengths and the complex conjugate branch does not ap- 
pear. For the shortest wavelengths and for all negative fir 
there are three real solutions: two negative solutions, one 
with the limit given by Eq. f75)) and the other approaching 
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Fig. 10. Real part of solutions of Eq. ([5^ as a function 
of the wave length l/{kH) for ao = 0.1, /3 = 0, ^ = 1 and 
various values of the time delay t. The bottom panel zooms 
iir the shaded area on the top parrel. The line converrtion 
is the same as iir previous figures. In this figure, we take 
7 = 5/3. 



— oo (Eq. [81]), arrd one positive solution approachiirg 
accordiirg to Eq. 



5. Summary 

We have derived the dispersion relation, Eq. ([5^ . for per- 
turbations of air accretion disk with heating that is offset 
in time relative to pressure perturbations. The standard a- 
prescriptioir was generalized to account for a time shift r be- 
tween the viscous stress response and perturbation of pres- 
sure, as well as for an arbitrary ratio of the corresponding 
perturbations fSect. l^TT]) . No restrictions were placed oir 
the allowed gas pressure to total pressure ratio, < /? < 1. 

In the limit of loirg waves the number of real solutioirs 
for the perturbatioir growth rate, fin, and their signs de- 
peird both on the relatioir between ^ and /3, aird on the 
value of the time lag — r. For all cases there is one triv- 
ial solution n = 0. For the standard case with no time lag 
(fir = 0) there is an additional real solution. It is negative 
if (Eq. EI]), 



The same condition applies wheir Qt > 0, i.e., advanced 
heating does irot affect the appearairce of the viscous in- 
stability in radiatioir pressure dominated disks. However, 
retarded heating may stabilize (or destabilize) the disk. For 
iir < 0, if inequality (|86|) is satisfied then two real roots 
with opposite signs exist (in additioir to the trivial one). If 
the inequality is irot satisfied two iregative roots appear, but 
only for delays larger than a critical value, Or < $7x2 < 0. 
The specific value of this critical flT2 is inversely propor- 
tional to ao and depends on f3 and ^ according to Eq. (|7D|) . 
For smaller values of the time lag, the non-trivial real solu- 
tions are positive or do not exist. Properties of real solutions 
of the dispersion relation in the limit of long wavelengths 
are summarized in Table [1] 

Iir addition to the real roots of the dispersion equation 
discussed above there exist (for fir ^ 0) an infinite num- 
ber of complex periodic solutions. Their real part, i.e., the 
growth rate, is opposite in sign to fir (for moderate values 
ofrjr). 

For very short waves there are two negative real so- 
lutions (one of the damping rates is finite, the other ap- 
proaches inifnity) independently of the system parameters. 
For rtr < there is an additional positive root diverging to 
-l-oo. 

Based on these properties we may conclude that the 
thermal and secular branches are stable for fir > if cri- 
terion is satisfied. For negative Ht the thermal branch 
is stable only if the time lag is sufRciently large. For neg- 
ative D,T complex solutions with a positive real part exist 
(Fig. [31) and therefore the issue of stability in this regime is 
more complicated. The growth rates of the secular branch 
are positive for all time lags (r < 0). 

In a forthcoming paper we shall present a detailed physi- 
cal discussion of the stability of disks with retarded heating, 
applying our conclusions to recent results of MHD numer- 
ical simulations of radiation-pressure dominated disks. 
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